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A  three-dimensional  and  two-phase  model  was  employed  to  investigate  the  effect  of  the  anisotropic 
GDL  thermal  conductivity  on  the  heat  transfer  and  liquid  water  removal  in  the  PEMFCs  with  serpentine 
flow  field  and  semi-counter  flow  operation.  The  GDL  with  different  anisotropic  thermal  conductivity  in 
the  three  directions  (x,  y,  z )  was  simulated  for  four  cases.  As  a  result,  the  water  saturation,  temperature, 
species,  current,  potential  distribution  and  proton  conductivity  were  obtained.  According  to  the  compar¬ 
ison  between  the  results  of  each  case,  some  new  conclusions  are  obtained  and  listed  as  below:  (1)  The 

anisotropic  GDL  produces  the  high  temperature  difference  than  that  of  isotropic  case,  and  the  in-plane 

thermal  conductivity  perpendicular  to  the  gas  channels  is  more  important  than  that  of  along  channels, 
which  may  produce  the  larger  temperature  difference.  (2)  Water  saturation  decreases  due  to  the  large 
temperature  difference  in  the  anisotropic  case,  but  some  water  vapor  may  condense  in  the  area  neighbor 
to  the  channel  ribs  due  to  the  cool  function  of  the  current  collector  and  the  great  temperature  difference. 

(3)  The  anisotropic  thermal  conductivity  in  the  through-plane  direction  and  the  in-plane  direction  per¬ 
pendicular  to  the  gas  channels  can  lead  to  the  decrease  of  the  membrane  conductivity.  (4)  The  isotropic 
GDL  is  better  than  that  of  anisotropic  one  for  the  uniform  current  density.  Also,  in-plane  thermal  con¬ 
ductivity  perpendicular  to  the  channels  has  more  negative  effect  on  the  current  density  distribution  in 
the  membrane  than  that  of  the  along  channels  one. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Heat  management  and  water  management  are  key  issues  in 
the  operation  of  proton  exchange  membrane  fuel  cells  (PEMFCs), 
which  are  well  recognized  by  the  researchers.  Also,  many  mod¬ 
els  and  experiments  have  been  done  to  investigate  and  optimize 
the  operation  to  achieve  the  better  water  management  and  good 
performance  of  PEMFCs  [1-6].  And  most  of  the  researches  are  con¬ 
centrated  on  the  effect  of  operating  conditions  (pressure,  flow  rate, 
stoichiometry,  humidified  level,  etc.)  on  the  water  removal,  how¬ 
ever,  the  heat  management  also  has  important  effect  on  the  water 
management  due  to  the  closely  interacted  thermal-mass  transfer 
coupled  behavior  between  water  and  temperature,  also  a  signif¬ 
icant  amount  of  heat  is  released  or  absorbed  due  to  the  phase 
change  of  water,  and  the  rates  of  which  are  also  a  strong  function  of 
temperature.  Therefore,  water  management  in  PEMFCs  should  be 
considered  with  thermal  management  simultaneously.  Hwang  [7] 
developed  a  two-phase  model  considering  velocity,  temperature 
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and  current  in  a  PEM  fuel  cell.  He  predicted  the  phase  equilib¬ 
rium  front  and  the  thermal  equilibrium  front  in  a  porous  cathode 
of  a  PEM  fuel  cell  and  concluded  that  increasing  the  rib-shoulder 
temperature  will  reduce  the  condensation  zones,  because  the  hot 
rib-shoulder  surfaces  increase  the  nearby  fluid-phase  temperature 
that  increases  the  saturation  pressure.  Meng  et  al.  [8  ]  also  presented 
a  non-isothermal  two-phase  model  of  PEMFCs,  and  their  results 
indicate  a  condensation/evaporation  interface  would  appear  in  the 
porous  materials  and  its  location  changes  with  the  inlet  humidity 
value  under  a  low-humidity  inlet  condition.  Also  liquid  water  is 
mainly  produced  in  the  GDL  in  two  regions;  one  is  near  the  current 
collecting  land  owing  to  the  low  temperature  and  another  further 
inside  the  GDL  but  still  away  from  the  catalyst  layer.  Falcao  et  al. 
[9]  recently  presented  a  simple,  using  low  CPU,  steady-state,  one¬ 
dimensional  model  accounting  for  coupled  heat  and  mass  transfer 
occurring  in  a  PEM  fuel  cell.  The  model  outputs  are  the  temper¬ 
ature  and  concentration  across  the  cell  and  the  water  content  in 
the  membrane.  The  model  predicts  reasonably  well  the  influence 
of  current  density  and  RH  on  the  net  water  transport  coefficient. 
Humidified  cathode  and  especially  humidified  anode  streams  are 
needed  to  avoid  the  membrane  dehydration,  particularly  at  high 
current  density.  The  above  researches  are  based  on  the  isotropic 
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Nomenclature 

a  water  activity 

C  molar  concentration  (mol  m-3 ) 

cp  specific  heat  capacity(J  kg-1  K-1 ) 

Cj  mass  transfer  coefficient  (s-1 ) 

D  diffusion  coefficient  (cm2  s-1 ) 

F  Faraday’s  constant 

hL  latent  heat  (kj  kg-1 ) 

/  current  density  magnitude  (A  nrr 2 ) 

tgef  reference  current  density  (A  m-2 ) 

J  reaction  rate 

I<  permeability  (cm2) 

k  thermal  conductivity  (W  m-1 1<-1 ) 

M  molar  mass  (kg  mol-1 ) 

n  number  of  electrons 

P  pressure  (Pa) 

R  universal  gas  constant  (J  mol-1 1<-1 ) 

ST  source  term  in  energy  equation 

Sw  source  term  in  liquid 

s  water  saturation 

T  temperature  (K) 

V  velocity  (cm  s-1) 

Voc  open  circuit  voltage  (V) 

y  mass  fraction 

Greek  letters 
0  potential  (V) 

s  volume  fraction 

a  surface  tension  (N  cm-1) 

<7mem  proton  conductivity  in  membrane 

crsoi  electric  conductivity  of  solid 

0C  equilibrium  contact  angle  on  diffuser 

A  polymer  water  content  H20/S03 

ad  water  drag  coefficient  in  membrane 

r]  overpotential  (V) 

fi  viscosity  (Pas) 

p  density  (kg cm-3) 

Subscripts 
an  anode 

C  about  capillary 

cat  cathode 

sat  Saturated 

sol  about  electron 

g  gas  phase 

H2  hydrogen 

i  note  for  species 

k  note  for  phases 

l  liquid  phase 

m  mixture  properties  of  multiphase  mixture 

mem  polymer  phase 

02  oxygen 

ohm  about  electronics 

p  phase 

q  phase 

ref  reference 

w  water 

wv  water  vapor 


GDL.  However  the  GDL  in  PEMFCs  made  of  carbon  paper  which 
has  anisotropic  properties  such  as  permeability,  thermal  conduc¬ 
tivity,  and  electronic  conductivity  [10].  So  the  research  based  on 
the  anisotropic  GDL  is  very  necessary. 


By  now,  also  some  researches  are  performed  to  investigate  the 
effect  of  the  anisotropic  properties  on  the  performance  of  PEM¬ 
FCs  [11-17].  The  influence  of  binder  structure  and  PTFE  treatment 
on  the  anisotropic,  effective  diffusivity  of  different  carbon  paper 
GDLs  has  been  experimentally  investigated  for  the  first  time  in 
Reto  Fliickiger’s  work  [11].  The  results  revealed  a  high  degree  of 
anisotropy  given  by  the  orientation  of  the  fibers  is  preferable  for 
homogeneous  conditions  under  flow  field  channel  and  rib.  It  is 
expected  to  reduce  degradation  effects  caused  by  local  current 
peaks.  Bapat  [12]  developed  a  two-dimensional  two-phase  model 
to  analyze  the  effects  of  anisotropic  electrical  resistivity  on  cur¬ 
rent  density  and  temperature  distribution  in  a  PEM  fuel  cell  and 
found  that  a  higher  in-plane  electrical  resistivity  of  the  gas  dif¬ 
fusion  layer  (GDL)  adversely  affects  the  current  density  in  the 
region  adjacent  to  the  gas  channel  and  generates  slightly  higher 
current  densities  in  the  region  adjacent  to  the  current  collector. 
He  also  developed  a  two-dimensional  two-phase  model  based  on 
the  classical  two-fluid  model  to  analyze  the  effect  of  low  through- 
plane  and  high  in-plane  thermal  conductivities  on  the  two-phase 
behavior  [13].  It  is  concluded  that  the  current  density  may  be  max¬ 
imized  at  low-humidity  operating  conditions  by  tailoring  the  GDL 
to  have  high  through-plane  thermal  conductivity  near  the  inlet 
and  progressively  decreasing  through-plane  thermal  conductivity 
at  distances  farther  away  from  the  inlet  along  the  flow  channel. 
Ramousse  [14]  developed  an  analytical  approach  to  estimate  the 
thermal  conductivity  of  GDL  composed  of  non-woven  carbon  felts 
with  anisotropicity.  Ju  [15]  presented  two-dimensional  simula¬ 
tion  to  study  the  effects  of  gas  diffusion  layer  (GDL)  anisotropy 
and  the  spatial  variation  of  contact  resistance  between  GDLs 
and  catalyst  layers  (CLs)  on  water  and  heat  transfer  in  polymer 
electrolyte  fuel  cells  (PEFCs).  The  simulation  results  clearly  demon¬ 
strate  that  GDL  anisotropy  and  the  spatial  variation  of  GDL/CL 
contact  resistance  have  a  strong  impact  on  thermal  and  two-phase 
transport  characteristics  in  a  PEFC  by  significantly  altering  the  tem¬ 
perature,  water  and  membrane  current  density  distributions,  as 
well  as  overall  cell  performance.  Hao  [16]  developed  a  multiple- 
relaxation-time  (MRT)  lattice  Boltzmann  method  (LBM)  to  predict 
the  anisotropic  permeability  of  GDL  with  consideration  of  porosity 
and  tortuosity.  In  Yang’s  work  [17],  a  coupled  electron  and  two- 
phase  mass  transport  model  for  anisotropic  GDLs  is  developed.  The 
effects  of  anisotropic  GDL  transport  properties  due  to  the  inherent 
anisotropic  carbon  fibers  and  caused  by  GDL  deformations  are  stud¬ 
ied.  Results  indicate  that  the  inherent  structural  anisotropy  of  the 
GDL  significantly  influences  the  local  distribution  of  both  cathode 
potential  and  current  density.  But  the  heat  transfer  is  not  included 
in  their  simulation. 

Although,  some  works  have  been  done  to  consider  the 
anisotropicity  of  the  GDL,  most  of  them  are  concentrated  on  the 
anisotropic  permeability,  however,  the  through-plane  permeabil¬ 
ity  is  not  greatly  different  from  the  in-plane  permeability  [10].  But, 
there  is  great  difference  between  through-plane  thermal  conduc¬ 
tivity  and  in-plane  thermal  conductivity  [14,15],  so  the  research 
of  the  anisotropic  thermal  conductivity  should  be  performed  due 
to  the  importance  of  the  heat  management,  and  only  a  few  refer¬ 
ences  are  reported  [13,1 5]  which  are  based  on  the  two-dimensional 
simulation,  and  obviously  different  from  the  actual  PEMFCs.  Fur¬ 
thermore,  the  effect  of  liquid  on  the  heat  transfer  is  not  considered 
which  can  be  confirmed  by  the  governing  equations  in  their  works 
[15].  Also,  the  effect  of  flow  field  on  the  heat  management,  and  the 
difference  in  the  two  directions  of  in-plane  thermal  conductivity 
cannot  be  involved  for  the  two-dimensional  cases. 

In  this  study,  a  two-phase  three-dimensional  simulation  is 
present  for  the  consideration  of  the  anisotropic  thermal  conduc¬ 
tivity  in  the  three  directions  (through-plane,  in-plane  along  the 
channel,  in-plane  perpendicular  to  the  channel)  of  the  GDLs  for  the 
commonly  used  serpentine  flow  field  with  semi-counter  flow  oper- 
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ation,  and  some  new  conclusions  are  obtained,  which  may  provide 
guidelines  for  the  design  and  operation  of  PEMFCs.  However,  gas 
diffusion  is  the  dominate  function  in  the  GDL,  and  the  convection 
is  very  weak  for  the  serpentine  flow  field,  so  the  effect  of  the  mass 
transfer  on  the  heat  transfer  is  limited  and  the  anisotropic  property 
for  mass  transfer  was  not  considered. 

2.  Three-dimensional  two-phase  simulation 

The  gas  is  treated  as  the  ideal  gas.  Due  to  the  relatively  low  veloc¬ 
ity,  the  fluid  flow  is  considered  as  laminar  flow.  The  liquid  contact 
angle  of  the  catalyst  layer  is  70°.  And  the  effect  of  the  channel  wall 
on  the  liquid  transfer  is  neglected. 

2.1.  Governing  equations 


Table  1 

Sources  terms. 

Source  terms  (zero  in  other 
region) 


Water  produced  (cathode 
catalyst  layer) 

Darcy  pressure  drop  of  gas 
in  cathode  and  anode, 

5  Drag 

Oxygen  reaction  rate 
(cathode  catalyst  layer) 
Hydrogen  reaction  rate 
(anode  catalyst  layer) 
Mass  transfer  rate  between 
gas  and  liquid 
Energy  source  term 


Defining  equation 


Sw=7c^ 

-£(i^?  +  ^)^ 


So2=-Jc-^ 
Sh2  =  ~Ja  -jjr 


r„  =  ct max  (  [(1  -  ,  [-sp,]) 

St  =  I2Rohm  +  Van,  cat  Jan,  cat  + 


The  continuity  equation  for  the  multiphase  mixture  is 

9 

+  v  ■  (pmVm)  =  0 

EL (skPkrt) 

Vm  —  - 

n 

Pm  =  TVi  {Skpk) 


where 

(1) 

■W 

ll 

hn 

(2) 

k= l 

keff  =  s , 

skPlJk 

Pm 


k=  1 


(12) 


The  momentum  equation  for  the  mixture  can  be  expressed  as  [18] 

9 

g^CPm^m)  +  V  -  (Pm^m^m) 

=  -Vp  +  V  •  [pm(V VmT  +  Vt/m)]  +  Pmg  T  Sdraag 

En 

,  [SkP-k)  (3) 

k= 1 


Pm 6pm  —  ^  ^Pk^p, 


k= 1 

For  the  anisotropic  GDL,  the  thermal  conductivity  can  be  described 
as 

(xx  0  0  ] 

=  0  yy  0  l  (13) 

0  0  zz 


From  the  continuity  equation  for  the  liquid  phase  water,  the  volume 
fraction  equation  for  liquid  water  can  be  obtained: 

9 

^(spi)  +  V-(esp,v£)  =  rw+Sw  (4) 


In  the  gas  diffusion  layer  and  catalyst  layer,  the  liquid  is  driven  by 
capillary  force,  then  it  can  be  rewritten  as 


V  ■  (sspivZ)  =  V  ■ 


—  rw  +  5y i 


(5) 


Capillary  pressure  is 

/  £  \  V2 

Pc=Pg-Pi=°  cosOl-j  J(s)  (6) 

The  species  mass  conservation  in  gas  phase  is 


v  ■  ( SgPgVgyt )  =  V  •  (DjVy,)  +  Si 


(7) 


The  membrane  phase  and  solid  phase  potential  conservation  equa¬ 
tions  and  electrochemical  reaction  rate  in  the  cathode  side  and 
anode  side: 


(  .  \  f  Ja  anode  catalyst  layer 

(crmem  Vcpmem )  -  j  cathode  catalyst  layer 

(8) 

,  ,  .  j  -Ja  anode  catalyst  layer 

asoi  ^so1  ~]^-Jc  cathode  catalyst  layer 

O) 

Water  motion  through  the  membrane  is 

V  -  (-DWVCW  +  A^)  =0 

(10) 

The  heat  transfer  equation  is 

9 

■^:(PmCp,m)  +  V  •  (PmCp,mVmT)  =  V  •  (/<ejjVT)  +  Sj 

(11) 

xx,  yy,  and  zz  are  the  direction  vectors. 

The  source  terms  and  parameters  in  the  equations  are  listed  in 
Tables  1  and  2. 

2.2.  Boundary  conditions  and  parameters 

The  simulation  domain  is  shown  in  Fig.  1,  which  includes  the 
current  collector,  gas  channel,  gas  diffusion  layer,  catalyst  layer 
both  in  anode  side  and  cathode  side,  and  membrane,  also  the  inlet 
and  outlet  for  the  anode  and  cathode  are  schematically  shown  in 
Fig.  1.  The  conventional  serpentine  flow  field  is  investigated  in  the 
present  study  (shown  in  Fig.  2),  which  has  1 7  channels.  And  the  size 
of  the  PEMFC  is  5  cm  x  3.3  cm. 

The  inlet  volume  flow  rate  is  500  ml  min-1,  which  is  converted 
to  mass  flow  rate  by  the  UDF  (user  defined  function  in  Fluent® 
software),  the  outlet  boundary  condition  is  the  pressure  outlet  con¬ 
dition,  the  outlet  pressure  is  equal  to  atmosphere  pressure.  The 
operation  current  density  is  1.0  A  cm-2.  The  inlet  temperature  is 


Anode  outlet 


Cathode  current 
collector 

CAthode  channel 
-MEA 

Anode  channel 
Anode  current  collector 


Cathode  outlet 


Fig.  1.  Calculation  domain. 
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Table  2 

Parameters. 

Parameters 

Reaction  rate  in  cathode  CL 

Reaction  rate  in  anode  CL 

Polymer  phase  conductivity 
Water  content  in  polymer  phase 
Water  index  in  polymer 
Water  activity 


Water  diffusivity  in  polymer  phase 


The  effective  diffusion  coefficient  [19] 

Leverett  J-function 

Saturated  water  vapor  pressure 


Defining  equation 


Jc  =  (l-*yW£^j  [exp(  = 

Ja  —  ( 1  —  sMvI'o,a  ^  ^  [exP  (  sol  ~  0mem  —  exP  ^  (0soJ  —  0mem  j 

&mem  =  £mem(5.1 4A  —  3.26)  exp  |^1 268  ^ —  y=^j(A,>l) 


r  -  _e+_ 

-  A+1 

A  =  0.043  +  17.81a  -  39.85a2  +  36.0a3 
k  =  14  +  1.4(a  -  1) 
a=pr+2  s 


a  <  1 
a  >  1 


Dw  =  10_10exp 
Dw  -  10“10exp 
Dw  =  10_10exp 
Dw  =  10“10exp 


2416  ( 

(  1 

-1)1 

^303 

T)  J 

2416  ( 

(  1 

-1) 

^303 

T)  \ 

2416  ( 

(  1 

-1) 

^303 

T)  J 

2416  ( 

(  1 

-1) 

^303 

T)  \ 

(2.563  -  0.33 k  +  0.0264A2  -  0.00067U3) 
(-l.25A.-h6.65) 

(2.05A  -  3.25) 


A.  >4 

3  <  A  <4 

2  <  A  <  3 

A  <2 


D|  =  £1,5(1  -s)2'5D,,o  (y )  (^)15 

f  1.417(1  -s)-  2.120(1  -s)2  +  1.263(1  -s)3  if0c  <  90° 

J{S  ~  \  1.417s- 2.120s2  +  1.263s3  if0c  >  90° 

log*"  =  -2.1794 +  0.02953(7-273.17) -9.1837  x  10“5(T  -  273. 17)2  +  1.4454  x  10“7(T  -  273.17)3 


353  K,  and  the  temperature  of  anode  current  collector  and  cathode 
current  collector  end  walls  boundary  are  the  353  K,  which  means 
the  current  collector  is  ideally  cooled  to  keep  constant  temperature. 
The  other  lateral  walls  and  the  end  walls  are  impermeable  for  all 
the  species.  The  potential  on  the  anode  current  collector  boundary 
is  set  to  be  0.  The  water  saturation  in  cathode  and  anode  inlet  are 
zero.  The  inlet  gas  for  anode  side  and  cathode  side  are  humidified 
at  348  K  before  entering  the  PEMFCs. 

The  values  of  the  other  parameters  used  in  the  model  are  listed 
in  the  Table  3. 

2.3.  Mesh  grid  and  solution  technique 

The  geometry  model  shown  in  Fig.  1  was  discretized  into 
660,000  hexahedral  mesh  volumes,  and  to  assure  the  quality  of 
the  grid,  the  size  of  the  grid  in  gas  channel,  gas  diffusion  layer, 
catalyst  layer  and  membrane  are  different.  The  simple  algorithm 
is  applied  for  solving  the  pressure-velocity  coupled  equations,  and 
species  equations.  And  suitable  relax  factors  are  used  for  water  sat¬ 
uration,  potential  and  water  content.  The  simulation  is  performed 
in  Flunet®6.3  software  of  Ansys  company  with  the  additional 
UDFs(User  Defined  Function  in  Fluent)  developed  by  the  authors. 

3.  Results  and  discussion 

To  investigate  the  effect  of  the  anisotropic  GDL  on  the  heat  and 
water  management,  the  simulations  for  four  cases  are  completed. 
The  conditions  of  the  four  cases  are  listed  in  Table  4.  And  the  typical 
anisotropic  values  of  the  thermal  conductivity  in  GDL  are  obtained 
in  Refs.  [20,14,15]. 

The  through-plane  thermal  conductivity  of  GDL  is  obtained  in 
Refs.  [20,14],  which  is  great  smaller  than  that  of  the  in-plane  case, 
this  is  determined  by  the  structure  of  the  carbon  paper.  In  com¬ 
mon  case,  the  carbon  paper  is  composed  of  random  packed  carbon 
fibers,  and  the  heat  transfer  in  the  through-plane  mainly  occurs 
between  the  different  carbon  fibers,  furthermore,  there  is  also  space 
between  the  carbon  fibers,  so  the  thermal  conductivity  is  weak.  As 
for  the  heat  transfer  in  the  in-plane  direction,  the  heat  transfer  is 
mainly  accomplished  in  the  same  carbon  fiber,  so  the  thermal  con¬ 


ductivity  is  larger  than  that  of  through-plane  case,  and  that  means 
the  in-plane  thermal  conductivity  is  determined  by  the  orientation 
and  length  of  the  carbon  fibers.  So,  the  thermal  conductivity  in  the 
in-plane  can  also  be  different.  In  the  present  study,  case  3  and  case 
4  are  designed  to  consider  this  point.  In  case  3,  it  means  that  all  of 
the  carbon  fibers  are  orientated  to  the  in-plane  direction  perpen¬ 
dicular  to  the  gas  channel,  and  case  4  means  all  of  the  carbon  fibers 
are  orientated  to  the  in-plane  direction  along  the  gas  channel.  So, 
the  thermal  conductivity  in  the  other  directions  is  the  same  as  the 
through-plane  one.  By  this  way,  the  effect  of  anisotropic  thermal 
conductivity  in  GDL  on  the  heat  transfer  is  related  to  the  structure 
of  the  flow  field. 

3.1.  The  effect  of  anisotropic  thermal  conductivity  on  the 
temperature 


Figs.  3-6  show  the  temperature  distribution  on  the  cathode 
catalyst  layer-membrane  surface  for  the  four  different  cases.  Corn- 


Fig.  2.  Schematic  view  of  the  flow  field. 


Table  3 

Values  of  the  parameters. 
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Physical  properties 

Value 

Faraday’s  constant,  F 

96487  C  mol-1 

Permeability  of  gas  diffusion  layer,  I< 

8  x  10-8  cm2 

Liquid  water  viscosity, 

3.565  xl0“4  Pas 

Anodic/cathodic  transfer  coefficient, aaloic 

0.5/0.55 

Water  contact  angle  in  diffuser,  0 

120° 

Gas  channel  width 

0.1  cm 

Gas  channel  height 

0.1  cm 

Thickness  of  current  collector 

0.15  cm 

Anode  GDL  thickness 

0.019  cm 

Cathode  GDL  thickness 

0.019  cm 

Gas  diffusion  layer  void  fraction 

0.78 

Catalyst  layer  thickness 

0.002  cm 

Catalyst  layer  void  fraction 

0.5 

Membrane  thickness  (Nation®  112) 

0.0005  cm 

Cell  Inlet  temperature 

353  K 

Outlet  pressure 

0.1  MPa 

Air  and  fuel  inlet  humidified  temperature 

348 1< 

Open  circuit  voltage 

0.95  V 

Mass  transfer  rate  between  phases,  cj 

100  s-1 

Gas  constant,  R 

8314Jkmol-1  K-1 

Reference  hydrogen  concentration,  CjJ 

1  kmol  m3 

Reference  oxygen  concentration,  cj 

1  kmol  m3 

Operation  current  density 

1.0  A  cm-2 

Anode  exchange  current  density, 

1.5e8  Am-3 

Cathode  exchange  current  density, 

7000  Am-3 

Thermal  conductivity  of  catalyst  layer  [20] 

0.27  Wm-1  K"1 

Thermal  conductivity  of  membrane  [20] 

0.29  Wm-1  K"1 

Specific  heat  capacity  of  membrane 

800 J  kg-1  K"1 

Specific  heat  capacity  of  GDL 

lOOOJkg-1  K"1 

Specific  heat  capacity  of  catalyst  layer 

lOOOJkg-1  K"1 

Specific  heat  capacity  of  current  collector 

800 J  kg-1  K"1 

Specific  heat  capacity  of  air 

1006.43  J  kg”1  K1 

Specific  heat  capacity  of  hydrogen 

14283  J  kg-1  K1 
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Fig.  3.  Temperature  on  the  cathode  catalyst  layer- membrane  interface  for  case  1, 
I<  ( x\m ,  z:m). 


Fig.  4.  Temperature  on  the  cathode  catalyst  layer-membrane  interface  for  case  2, 
K  ( x:m ,  z:m). 
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Fig.  5.  Temperature  on  the  cathode  catalyst  layer-membrane  interface  for  case  3, 
K  ( x:m ,  z:m). 


pared  with  the  temperature  for  the  isotropic  case,  it  is  obvious  that 
the  maximum  temperature  for  all  the  anisotropic  cases  (357.6  °C, 
360.4  °C,  and  357.8  °C  separately)  are  larger  than  the  maximum 
temperature  for  the  isotropic  case  (356.4  °C).  This  is  due  to  the  weak 
thermal  conductivity  of  the  anisotropic  cases  (the  thermal  conduc¬ 
tivity  for  the  isotropic  case  in  the  through-plane  is  nearly  10  times 
of  that  for  the  anisotropic  cases).  The  maximum  temperatures  are 
also  different  for  the  different  anisotropic  cases.  And  according  to 
the  result  in  Figs.  4-6,  there  is  great  difference  for  the  temper¬ 
ature  distribution  for  the  different  anisotropic  cases.  The  largest 


Table  4 

Different  thermal  conductivity  for  the  simulation. 


Through-plane  thermal 
conductivity  (y)  (W  m_1  K1 ) 

In-plane  thermal  conductivity 
along  channels  (z)  (W  m_1  K-1 ) 

In-plane  thermal  conductivity  perpendicular  to 
the  channels  (x)  (W  irr1  K_1 ) 

Case  1 

10  [15] 

10 

10 

Case  2 

1.27  [14,20] 

10 

10 

Case  3 

1.27 

10 

1.27 

Case  4 

1.27 

1.27 

10 
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Fig.  6.  Temperature  on  the  cathode  catalyst  layer-membrane  interface  for  case  4, 
I<  (x:m,  z:m). 


temperature  difference  between  the  area  neighbor  to  the  channel 
area  and  the  area  neighbor  to  the  ribs  is  obtained  in  case  3.  This  is 
caused  by  the  heat  in  the  area  neighbor  to  the  channel  area  is  hard 
to  transfer  due  to  the  weak  thermal  conductivity  in  the  x  direc¬ 
tion,  which  makes  the  temperature  increase  sharply.  However,  For 
case  2  and  case  4,  the  temperature  distribution  including  the  max¬ 
imum  temperature  are  nearly  the  same  (357.6  °C  and  357.8  °C  for 
the  maximum  temperature).  That  means  the  weak  in-plane  along- 
channel  thermal  conductivity  has  little  effect  on  the  heat  transfer  in 
the  fuel  cells.  However,  all  the  operation  conditions  and  the  prop¬ 
erties  of  the  material  are  the  same  for  case  3  and  case  4,  and  the 
only  difference  is  the  direction  of  the  weak  thermal  conductivity, 
for  case  3,  through  plane  of  and  in-plane  perpendicular  to  the  gas 
channel  thermal  conductivities  are  weak,  while  the  weak  thermal 
conductivity  are  in  through  plane  of  GDL  and  in-plane  along  the  gas 
channel  (except  the  corner)  for  case  4.  In  the  present  fuel  cell,  the 
heat  produced  in  the  operation  is  removed  in  the  following  ways: 
firstly,  the  heat  is  removed  through  the  anode  and  cathode  current 
collector  which  neighbors  the  cool  plate  in  the  actual  application, 
the  function  of  the  cool  plate  is  described  by  the  constant  tem¬ 
perature  on  the  current  collector-cool  plate  interface  (end  walls  of 
current  collector)  in  the  present  study.  In  this  case,  the  heat  trans¬ 
fer  mainly  happened  in  the  direction  of  through  plane,  but  there 
is  no  difference  between  case  3  and  case  4  related  to  that,  so  the 
heat  transfer  in  this  way  should  also  be  the  same.  The  second  way 
of  heat  removal  is  heat  removed  by  the  reactant  gas,  which  flows 
in  fuel  cell  and  flows  out  carrying  the  excessive  heat.  Because  of 
the  higher  thermal  conductivity  of  the  current  collector  compared 
to  that  of  the  reacting  gases,  heat  can  be  easily  removed  from  the 
current  collector,  so  the  temperature  in  the  area  near  the  ribs  of  the 
current  collector  is  low  compared  to  that  of  the  area  near  the  chan¬ 
nel,  which  can  be  demonstrated  by  all  the  simulation  temperature 
distribution.  For  case  3,  the  in-plane  thermal  conductivity  perpen¬ 
dicular  to  the  gas  channel  is  weak,  that  means  the  heat  transfer  from 
the  GDL  area  near  gas  channels  to  the  GDL  area  near  ribs  is  hard, 
and  the  heat  transfer  ability  in  the  area  near  gas  channels  is  also  low 
due  to  the  weak  thermal  conductivity  and  the  limit  flow  rate  of  the 
gases.  So,  the  produced  heat  cannot  be  removed  efficiently,  which 
causes  the  high  temperature  near  the  channels  and  the  great  tem¬ 
perature  difference  between  area  near  channels  and  area  near  ribs. 
In  other  words,  in-plane  thermal  conductivity  perpendicular  to  the 
gas  channels  is  very  important  for  the  heat  transfer  in  the  PEMFCs. 


As  for  the  case  4,  the  weak  thermal  conductivities  are  through- 
plane  and  in-plane  along-channel  thermal  conductivities.  But  there 
is  little  difference  between  the  results  of  case  2  and  case  4.  That 
means  the  heat  transfer  is  mainly  accomplished  in  the  through- 
plane  direction  and  in-plane  direction  perpendicular  to  the  gas 
channels,  so  the  in-plane  thermal  conductivity  along  the  gas  chan¬ 
nels  is  not  important,  and  has  limited  effect  on  the  heat  transfer  in 
the  PEMFCs.  So,  it  can  be  concluded  that  the  anisotropic  GDL  pro¬ 
duces  the  high  temperature  difference  than  that  of  isotropic  case, 
and  the  in-plane  conductivity  perpendicular  to  the  gas  channels  is 
more  important  than  that  of  along  channels,  and  the  anisotropic 
thermal  conductivity  in  this  direction  may  produce  the  even  larger 
temperature  difference. 


3.2.  The  effect  of  anisotropic  thermal  conductivity  on  the  water 
saturation  in  cathode  GDL  and  catalyst  layer 

Figs.  7-10  show  the  water  saturation  in  cathode  catalyst  layer 
and  gas  diffusion  layer  in  the  cross-section  (z  =  0.025  cm)  of  the 
PEMFC.  It  can  be  seen  that  the  maximum  water  saturation  for  the 
isotropic  case  is  about  0.07,  while  the  maximum  value  for  the  other 
cases  are  about  0.05.  According  to  the  analysis  in  the  Section  3.1 
that  the  isotropic  GDL  has  the  low  temperature,  that  means  low 
saturated  water  vapor  concentration,  so  less  water  produced  by 
the  electrochemical  reaction  can  vaporize  to  the  water  vapor  com¬ 
pared  with  the  other  cases,  and  more  water  remains  in  liquid  which 
leads  to  the  high  water  saturation.  In  other  words,  the  anisotropic 
GDL  decreases  the  water  saturation  in  the  cathode  catalyst  layer 
and  GDL. 

It  is  well  known  that  the  liquid  water  flows  through  the  catalyst 
layer  and  GDL  is  driven  by  the  capillary  force,  and  there  is  water 
saturation  gradient  through  the  thickness  of  the  catalyst  layer  and 
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Fig.  7.  Water  saturation  in  the  cross-section  of  z  =  0.025m  of  cathode  GDL  and  cat¬ 
alyst  layer  for  case  1  ( x:m ,  z:m). 
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Fig.  8.  Water  saturation  in  the  cross-section  of  z  =  0.025m  of  cathode  GDL  and  cat¬ 
alyst  layer  for  case  2  ( x:m ,  z:m). 


the  GDL,  i.e.  water  saturation  increases  from  GDL-channel  surface 
to  the  catalyst  layer  surface,  which  is  confirmed  by  many  simula¬ 
tion  results.  And,  the  present  simulation  results  for  the  isotropic 
GDL  also  show  the  same  trend  with  described  above  which  can  be 
seen  in  Fig.  7.  But,  as  for  the  result  for  the  other  three  cases,  it  is 
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Fig.  9.  Water  saturation  in  the  cross-section  of  z  =  0.025m  of  cathode  GDL  and  cat¬ 
alyst  layer  for  case  3  ( x:m ,  z:m). 
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Fig.  10.  Water  saturation  in  the  cross-section  of  z  =  0.025m  of  cathode  GDL  and 
catalyst  layer  for  case  4  ( x:m ,  z:m). 


different.  It  can  be  seen  that  there  are  relative  high  water  satura¬ 
tion  areas  neighbor  to  the  current  collector’s  ribs,  not  only  in  the 
catalyst  layer  in  Figs.  8  and  10  for  case  2  and  case  4,  but  also  more 
obvious  in  Fig.  9  for  case  3,  that  means  water  saturation  neighbor 
to  the  ribs  is  larger  than  that  of  the  area  neighbor  to  the  catalyst 
layer,  and  the  motion  of  the  water  is  not  mainly  determined  by  the 
capillary  force,  there  is  also  something  else:  temperature.  Accord¬ 
ing  to  the  results  in  Figs.  3-6,  the  temperature  in  the  catalyst  layer 
and  GDL  of  case  2,  case  3  and  case  4  are  higher  than  that  of  case 
1  both  for  the  area  neighbor  to  the  ribs  and  neighbor  to  the  chan¬ 
nel.  So,  more  liquid  water  vaporizes  to  water  vapor  in  those  three 
cases  in  the  catalyst  layer  near  the  channels,  and  the  water  vapor  in 
the  area  neighbor  to  the  gas  channel  can  directly  enter  into  the  gas 
channel  and  flows  out,  which  makes  the  water  saturation  decrease 
from  the  catalyst  layer  to  the  gas  channel,  as  for  the  water  in  the 
area  neighbor  to  the  ribs,  also  more  water  vaporizes  to  water  vapor 
in  the  catalyst  layer  due  to  the  relative  high  temperature,  but  there 
is  long  distance  from  this  area  to  the  gas  channels,  which  leads  to 
the  high  flow  resistance  compared  with  the  former  one,  further¬ 
more  the  temperature  in  the  area  neighbor  to  the  ribs  is  low  due  to 
the  good  heat  removal  performance  of  the  current  collector,  which 
makes  some  water  vapor  from  the  catalyst  layer  condense  in  the 
area  neighbor  to  the  ribs,  so  the  water  saturation  increases  from 
catalyst  layer  to  the  GDL.  So,  it  can  be  concluded  that  water  sat¬ 
uration  decreases  due  to  the  large  temperature  difference  in  the 
anisotropic  case,  and  some  water  vaporized  in  the  catalyst  layer 
because  of  the  high  temperature  for  anisotropic  cases,  but  some 
water  vapor  condenses  in  the  area  neighbor  to  the  ribs  due  to  the 
well  cool  function  of  the  current  collector,  which  cause  the  water 
saturation  in  the  area  neighbor  to  the  ribs  is  larger  than  that  of  in 
the  catalyst  layer. 

3.3.  The  effect  of  anisotropic  thermal  conductivity  on  the  proton 
conductivity 

The  proton  conductivity  of  the  membrane  is  determined  by  the 
water  content  and  temperature  in  the  membrane,  high  tempera¬ 
ture  and  high  water  content  are  helpful  to  the  proton  conductivity. 
According  to  the  analysis  before,  the  temperature  increases  due  to 
the  anisotropic  thermal  conductivity  in  the  GDL,  but  as  the  tern- 
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Fig.  11.  Proton  conductivity  on  the  cathode  side  of  the  membrane  for  case  1,  S  m-1 
(x:m,  z:m). 


Fig.  13.  Proton  conductivity  on  the  cathode  side  of  the  membrane  for  case  3,  S  m-1 
( x:m ,  z:m). 


perature  increases,  some  liquid  water  vaporized  to  water  vapor, 
so  the  water  saturation  for  the  anisotropic  case  decreases  accord¬ 
ing  to  the  conclusion  in  Section  3.2.  To  investigate  the  effect  of  the 
anisotropic  thermal  conductivity  on  the  performance  of  the  mem¬ 
brane,  the  proton  conductivity  on  the  membrane-cathode  surface 
for  the  different  cases  is  shown  in  Figs.  11-14.  It  can  be  easily  seen 
that  the  proton  conductivity  of  the  case  3  (maximum  value  is  12.0, 
and  the  minimum  value  is  4.5)  is  smaller  than  that  of  case  2  and 
case  4  (maximum  value  is  12.0,  and  the  minimum  value  is  5.0). 
As  for  the  case  1  (the  isotropic  case),  the  maximum  value  is  12.0, 
while  the  minimum  value  is  6.0.  Also,  the  average  proton  conduc¬ 
tivity  of  the  membrane  for  the  four  cases  is  obtained  and  shown  in 
Fig.  15.  It  can  be  easily  seen  that  the  highest  proton  conductivity  is 
in  case  1  and  lowest  proton  conductivity  is  in  case  3,  while  the  pro¬ 
ton  conductivity  for  case  2  and  case  4  are  nearly  the  same,  which 
are  higher  than  that  of  case  3  and  lower  than  that  of  case  1 .  So,  the 
best  membrane  performance  is  obtained  for  isotropic  GDL,  and  the 
conclusion  can  be  drawn  as:  the  anisotropic  thermal  conductivity  in 
the  through-plane  direction  and  the  in-plane  direction  perpendic- 
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Fig.  14.  Proton  conductivity  on  the  cathode  side  of  the  membrane  for  case  4,  S  irr 1 
( x:m ,  z:m). 
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Fig.  12.  Proton  conductivity  on  the  cathode  side  of  the  membrane  for  case  2,  S  m-1 
(x:m,  z:m). 


Fig.  15.  Comparison  between  the  average  proton  conductivity  of  the  membrane  for 
four  cases. 
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Fig.  16.  The  current  density  in  the  middle  membrane  for  the  case  1,  Am-2  ( x:m , 
z:m). 


Fig.  18.  The  current  density  in  the  middle  membrane  for  the  case  3,  Am-2  ( x:m , 
z:m). 


ular  to  the  gas  channels  can  lead  to  the  decrease  of  the  membrane 
performance. 

3.4.  The  effect  of  the  anisotropic  thermal  conductivity  on  the 
current  density  in  the  membrane 

Figs.  16-19  show  the  current  density  in  the  middle  membrane 
for  the  four  cases.  It  can  be  seen  that  the  range  of  current  density 
in  the  membrane  for  the  isotropic  case  is  smaller  than  that  of  other 
three  cases,  and  the  current  density  difference  between  the  area 
near  the  ribs  and  the  area  near  channels  for  the  isotropic  case  is 
also  smaller  than  that  of  the  other  three  cases.  The  current  density 
distribution  in  the  middle  area  of  first  case  in  Fig.  16  is  nearly  uni¬ 
form,  but  for  the  other  three  cases,  the  difference  is  very  obvious, 
especially  for  the  third  case.  And  the  distribution  for  the  case  2  and 
case  4  are  nearly  same.  In  the  operation  of  PEMFCs,  the  uniform  cur¬ 
rent  density  distribution  is  pursued  due  to  the  non-uniform  current 
density  may  cause  the  distortion  and  other  things.  So,  in  the  view  of 
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Fig.  17.  The  current  density  in  the  middle  membrane  for  the  case  2,  Am-2  (x:m, 
z:m). 


Fig.  19.  The  current  density  in  the  middle  membrane  for  the  case  4,  Am-2  (x:m, 
z\m). 

this  point,  the  isotropic  GDL  is  better  than  that  of  anisotropic  one. 
Also,  in-plane  thermal  conductivity  perpendicular  to  the  channels 
has  more  negative  effect  on  the  current  density  distribution  in  the 
membrane  than  that  of  the  along  channels  one. 

4.  Conclusions 

In  this  paper,  a  three-dimensional  and  two-phase  model  was 
employed  to  investigate  the  effect  of  the  anisotropic  GDL  ther¬ 
mal  conductivity  on  the  heat  transfer  and  liquid  water  removal 
in  the  PEMFCs  with  serpentine  flow  field  and  semi-counter  flow 
operation.  The  GDL  with  different  anisotropic  thermal  conductiv¬ 
ity  in  the  three  directions  (x,  y,  z)  was  simulated  in  four  cases 
i.e.  case  1  (isotropic  thermal  conductivity  with  the  good  thermal 
conductivity  in  every  direction),  case  2  (anisotropic  through-plane 
thermal  conductivity,  which  means  the  carbon  fibers  are  packed  in 
the  through-plane  direction,  so  the  thermal  conductivity  in  this 
direction  is  weak),  case  3  (anisotropic  thermal  conductivity  in 
through-plane  direction  and  in-plane  direction  perpendicular  to 
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the  gas  channel,  which  means  all  the  packed  carbon  fibers  are  ori¬ 
entated  to  the  along-channel  direction,  so  the  thermal  conductivity 
in  the  other  two  directions  are  weak),  case  4  (anisotropic  ther¬ 
mal  conductivity  in  through-plane  direction  and  in-plane  direction 
along  to  the  gas  channel,  which  means  all  the  packed  carbon  fibers 
are  orientated  to  the  direction  perpendicular  to  the  gas  chan¬ 
nel,  so  the  thermal  conductivity  in  the  other  two  directions  are 
weak).  As  a  result,  the  water  saturation,  temperature,  species,  cur¬ 
rent,  potential  distribution,  water  content,  proton  conductivity,  etc. 
were  obtained.  The  effect  of  the  anisotropic  thermal  conductivity 
was  investigated  by  the  comparison  between  the  results  of  the  tem¬ 
perature,  water  saturation,  proton  conductivity  and  outlet  water 
ratio  for  the  four  cases.  And  the  results  show  that: 

(1 )  The  anisotropic  GDL  produces  the  high  temperature  difference 
than  that  of  isotropic  case,  and  the  in-plane  thermal  conductiv¬ 
ity  perpendicular  to  the  gas  channels  is  more  important  than 
that  of  along  channels.  The  anisotropic  thermal  conductivity  in 
this  direction  may  produce  the  larger  temperature  difference. 

(2)  Water  saturation  decreases  due  to  the  large  temperature  differ¬ 
ence  in  the  anisotropic  case,  and  some  liquid  water  vaporizes 
in  the  catalyst  layer  because  of  the  high  temperature  for 
anisotropic  cases,  but  some  water  vapor  condenses  in  the  area 
neighbor  to  the  ribs  due  to  the  well  cool  function  of  the  current 
collector  especially  for  the  anisotropic  cases,  which  makes  the 
water  saturation  neighbor  to  the  ribs  is  larger  than  that  of  the 
catalyst  layer. 

(3)  The  anisotropic  thermal  conductivity  in  the  through-plane 
direction  and  the  in-plane  perpendicular  to  the  gas  channels 
direction  can  lead  to  the  decrease  of  the  membrane  perfor¬ 
mance. 


(4)  In  the  view  of  current  density  distribution  in  the  membrane, 
the  isotropic  GDL  is  better  than  that  of  anisotropic  one.  Also, 
in-plane  thermal  conductivity  perpendicular  to  the  channels 
has  more  negative  effect  on  the  current  density  distribution  in 
the  membrane  than  that  of  the  along  channels  one. 

These  conclusions  provide  the  new  guidelines  for  the  design  of 
GDL  and  operation  of  PEMFCs. 
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